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We demonstrate for the first time a tight binding model for water incorporating polarizable anions. A novel 
aspect is that we adopt a "ground up" approach in that properties of the monomer and dimer only are 
fitted. Subsequently we make predictions of the structure and properties of hexamer clusters, ice-XI and 
liquid water. A particular feature, missing in current tight binding and semiempirical hamiltonians, is that 
we reproduce the almost two-fold increase in molecular dipole moment as clusters are built up towards the 
limit of bulk liquid. We concentrate on properties of liquid water which are very well rendered in comparison 
with experiment and published density functional calculations. Finally we comment on the question of the 
contrasting densities of water and ice which is central to an understanding of the subtleties of the hydrogen 
bond. 



I. INTRODUCTION 

We are motivated to construct a model of water for use 
in atomistic simulation within the framework of the self 
consistent polarizable ion tight binding (TB) theory 1,2 
The intention is that this should bridge the outstanding 
gap between the local density functional theory (DFT) 
and classical molecular mechanics models. The risk is 
that the model will prove too slow computationally to 
replace existing, highly effective classical models while 
being unable to capture enough of the physics and chem- 
istry (especially of the hydrogen bond) to provide a useful 
imitation of the first principles approach. To address the 
latter concern we have adopted a "ground up" strategy in 
which parameters of the TB model are fitted first to the 
monomer and lastly to the dimer, so that all the proper- 
ties of the solid and liquid are predictions of the model. 
The former concern may be settled by observing at once 
that our model is capable of propagating 500 self consis- 
tent molecular dynamics time steps for a unit cell of 128 
water molecules in one CPU hour over four threads of a 
commodity multicore 2.4GHz processor — 12 ps per day 
using a time step of 1 fs. 

The structure of the paper is this. We begin in sec- 
tion II by describing TB models constructed intuitively 
using parameters that are selected with little adjust- 
ment from published TB models, in particular drawing 
on Walter Harrison's seminal scheme. 3 This should give 
the reader confidence that the TB approach contains the 
essential physics. In section III we demonstrate a model 
that is contructed by tuning the parameters using a com- 
plex genetic algorithm that provides us with the most 
accurate rendering of the properties of the monomer, in 
particular its polarisability and force constants, and the 
dimer, in particular its geometry and the shape of the 
energy versus length of the hydrogen bond as deduced 
from high level quantum chemical calculations. We ap- 
ply this model first to the structure and energetics of the 
water hexamer and to ice (sections IV and V); and in sec- 
tion VI we demonstrate its predictive power in respect of 
the properties of liquid water, especially the radial dis- 
tribution function, dielectric constant and self diffusion 
coefficient. We discuss our results in section VII in which 



we also consider the question of the relative densities of 
ice and water and we conclude in section VIII. 



II. INTUITIVE MODELS 

The self consistent polarizable ion tight binding the- 
ory was proposed in the first instance to address the 
atomic and electronic structure of ceramic oxides whose 
properties are dominated by the highly polarizable oxy- 
gen anion. 1,4 Later it was applied to problems involv- 
ing polarizability in molecular physics. 2 The method is 
described in detail in a recent textbook 5 and lecture 
article. 6 We should point out a close similarity between 
our theory and that of the SCC-DFTB method. 7 The 
principal difference is that we go beyond the point charge 
approximation and allow the self consistent development 
of point charge multipoles; this is accompanied by associ- 
ated higher angular momentum components of the elec- 
trostatic potential allowing us to describe crystal field 
effects in a self consistent manner. The second difference 
is in our approach to the choice of model parameters, 
which we do largely by rational intuition rather than di- 
rect computation and fitting to density functional total 
energies. In section III when we demonstrate our "ge- 
netic" model, we will emphasize that for the special case 
of water, it is better to fit to high level quantum chem- 
istry and experiment than to DFT. 

A. Monomer 

1. A simple non self consistent model 

One might construct a very simple, non self consistent 
model for the monomer in the spirit of Harrison's "solid 
state table". 3 In this way we place a Is orbital on the 
hydrogen atoms and 2s and 2p orbitals on the oxygen 
with off diagonal hamiltonian matrix elements, or hop- 
ping integrals being 

V U , x = a m , x —- 
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in which r) ssa = — 1.4 and rj spa = 1.84. Here, r is the 
bond length and m the electron mass. \ stands for cither 
n or 7r bonding. The factor a is ours, since we have 
have found that a = \ produces a HOMO-LUMO gap in 
better agreement with experiment. Specifically in order 
to supress intermolecular O-H interactions and generally 
to achieve short ranged hopping, we adopt a modified 
scaling known as GSP, after its authors, 8 such that 



V(r)=V(ro)(^) n exp 



(1) 

The virtue of the GSP scaling is that is has the same 
value at ro (but not slope) as the Harrison form if we 
use n — 2 while it decays exponentially with distance 
under the control of the cut off parameters, r c and n c ; 
r is conventionally the equilibrium bond length, but as 
with all the parameters this could depend on the orbital. 
For this reason there are implicit {££'x} subscripts on 
the parameters in (1) that we have suppressed. In this 
intuitive model we set n = 2, n c = 4 and r c = 1.8r ; 
ro is our target bond length, 1.8094 bohr. The diago- 
nal matrix elements of the hamiltonian in the solid state 
table 3 are taken to be Hartree-Fock term values, namely 
e s = — 1 Ry on the H-atom, the values for oxygen being 
e s = -2.142 Ry and e p = -1.038 Ry. We should point 
out that in this, and indeed all our subsequent models, 
the bond angle is dominated by the s — p splitting on 
the oxygen atom. This is clear since in the absence of 
on-site s — p hybridisation the bond angle would be ex- 
actly 90°. Then as e s is raised relative to e p the bond 
angle increases accordingly. 9 It is notable that the bond 
angle is very accurately rendered using the atomic term 
values. 



2. Intuitive point charge model 



FIG. 1. A cartoon of the water monomer. The target bond 
angle is 6 = 104.26°, and the bond length is ro = 1.809 bohr. 
In (a), the dipole moment arising from charge transfer is indi- 
cated and it is a simple matter to adjust the charge transfer, 
S, by a choice of Hubbard-(7 to achieve the target dipole 
moment, 1.86 D. If we admit that the anion is polarizable 
(b) then the self consistent problem can also be solved. The 
charge transfer is increased by reducing U so that a larger 
charge transfer derived dipole moment exists which is in turn 
partially canceled by the induced dipole, again to achieve the 
target dipole moment. See the text for details of how this is 
done. 



q = Se 





(b) 



and stick to a non self consistent model. 10 On the other 
hand the self consistent inclusion of charge transfer is 
exactly what is needed to describe mixed covalent ionic 
bonding as we expect to find in water. Therefore at the 
next level of an intuitive model we specify Hubbard-Z7 
parameters which we choose to be 1 Ry on the H-atoms 
and which we adjust on the oxygen until the self consis- 
tent dipole moment has the value, 1.86 Debye (D), that is 
observed. We can now summarise this, our point charge 
model, by reference to Fig. 1(a). A fractional amount, 
S of an electron is transferred from each H-atom to the 
anion; the resulting charges give rise to a dipole moment, 



We will find that charge transfer which leads to an 
electrostatic repulsion between the H-atoms has a very 
small influence on the bond angle. Turning to the charge 
transfer, this is where the self consistency appears at the 
simplest level. After solution of the Schrodinger equation 
using the Harrison TB hamiltonian, examination of the 
eigenvectors will reveal Mulliken charges at the atomic 
sites. One might ignore these, or compute their contribu- 
tion to the energy using elementary electrostatics. Since 
charge transfer will always serve to lower the energy, if 
the electrostatic potential arising due to Mulliken charges 
is added to the hamiltonian and a further solution to the 
Schrodinger equation is made then, as this process is re- 
peated, charge will continue to accumulate at the more 
electronegative site resulting in a Coulomb catastrophe. 
This is corrected by a term in the hamiltonian, an energy 
proportional to the charge squared, such that the inter- 
site Coulomb energy is balanced by the on-site Hubbard 
energy so called. 10 These two are expected largely to can- 
cel in many cases so that it is admissible to neglect these 



p(S) = 2r (5ecos±<9 (2) 

which, once we have fixed our target bond angle, 9, and 
bond length, depends only upon <5 as p(S) = 5.655 D. 
Hence since we know our target dipole moment it is an 
easy matter to adjust our one free parameter, the anion 
Hubbard-[/, to obtain a point charge TB model for the 
electronic structure. There remains to fix the standard 
pair potential, which in the language of the tight bind- 
ing bond model is said to account for the all terms ex- 
cept for the bandstructure energy in the Harris-Foulkcs 
functional. 6,11 We adopt the GSP form (1) for this al- 
lowing two free parameters, namely, the prefactor which 
we call A(ro ) and the exponent n, which in the case of a 
pair potential we call m. We fit these exactly to the bond 
length and symmetric vibrational force constant from ex- 
periment, while setting m c = 6 and r c = 2.9 bohr. This 
completes the description of our point charge model. 



3 



3. Intuitive dipole model 

Let us now admit that the anion is a polarizable species 
and the dipole moment created by the charge transfer will 
result in an electric field acting at the oxygen atomic site. 
This will induce a dipole moment which we model as a 
point dipole 

2eS 

Pind = -ao (S) — cos ±0 (3) 
r o 

in which ao is a notional polarizability of the oxygen an- 
ion. If we knew what this was, and in principle it should 
depend on the oxidation state of the anion (as indicated 
by its dependence upon S), then our dipole model would 
follow immediately without any further guessing of pa- 
rameters. We would write using (2) and (3) 

Pind _ ct (5) 
p(S) ~ rg 

while 

P(<5)+Pind = l-86D 

which could be solved if we knew ao (<*>)■ In fact exam- 
ining the known polarizabilities of the neutral first row 
atoms, nitrogen through neon, we found it to be a rather 
linear function of the valence leading to an interpolation 
formula, 

a (S) w 0.8-0.46(5 (A 3 ) (4) 

Unfortunately using the resulting value of ao leads to a 
polarization catastrophe: The bond angle tends to zero 
to allow the electric field to maximize and induce the 
largest possible dipole moment at the anion. Therefore, 
to construct a sensible intuitive model, we need to guess 

a smaller polarizability and we choose 0.25 A There- 
after it should be straight forward to readjust the anion 
Hubbard-C/ until the target dipole moment is achieved. 
However the bond angle is reduced in this process from 
the same origin that gives rise to the polarization catas- 
trophe and it is necessary to raise the oxygen e s . But it is 
simple enough to adjust these two parameters simultane- 
ously, after which fitting the pair potential to r and the 
symmetric vibrational force constant finalizes the form of 
our intuitive dipole model. 

4. Crystal field 

We should make it clear at this point that our self 
consistent TB theory does not deal directly with anion 
polarizability. Instead the action of the electric field at 
an atomic site is captured by a term in the hamiltonian 
in addition to the non self consistent hamiltonian such 
that we have 

H = H + H' 



and H' amounts to a self consistent field having the form 

HrL' RL" = Ur <7r Sl'L" + ^ Vrl ^fJfJ'Z Cl'L"L (5) 

L 

and the polarizability is expressed through new parame- 
ters, A.£>(n( which are generalizations of the parameters 
(r l ) of crystal field theory. 12,13 In (5) R labels atomic 
sites and L is a composite angular momentum index 
L = {£m}; qr is the Mulliken charge transfer at site 
R and Ur is the Hubbard-i7 at that site. The second 
term acts to split the otherwise degenerate atomic energy 
levels £t through the appearance of on-site, off diagonal 
matrix elements of the hamiltonian. Cl'L"L arc Gaunt 
integrals (Al) and Vrl are expansion coefficients of the 
electrostatic potential (energy) seen by an electron at site 
R, this potential having been expanded into spherical 
waves 14 

Vr(v)=J2Vrl r e Y L (r) 

L 

Poisson's equation relates the electrostatic potential to 
the multipole moments Qr>l> of the Mulliken charge 
transfer 

Vrl = e 2 ^ Brlr'L' Qr'L' 

R/L' 

in which B is a sort of generalized Madelung matrix, 
equal to |R — R'| _1 for point charges. The general case 
is described in Appendix A. Finally, the multipole mo- 
ments are obtained from the self consistent eigenvectors, 
c, of the hamiltionian 1 ' 5 ' 6 

QRL = ^ X! f n CRL> c RL"Ae>l"l C L 'L"L 
L'L" n 

where /„ is the occupation number of state n. It is no- 
table that the crystal field parameter plays a dual role as 
multipole strengths in the theory. 

For the case of our water model the only polarizability 
parameter we need is the quantity A spp at the anion site. 
The interpretation of &inm is that the ^-component of 
the potential (in the case 1=1 the electric field) causes 
an on-site coupling in the hamiltonian of the i 1 and t" 
orbitals (in this case the s and p are coupled by the crys- 
tal field). In its role as multipole strength parameter, it 
describes the contributions of the on-site i 1 and i" com- 
ponents of the eigenvectors to the development of the 
£-pole moment of the charge (in this case the dipole). 
This reflects the well-known requirement to mix s and p 
orbitals on-site if one is to develop a dipole moment of 
the charge. 

5. Predicted properties of the monomer 

The question now is, what is the predictive power 
of these two models? Table I gives an answer. For 
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TABLE I. Predictions of the intuitive point charge and dipole 
models for the water monomer. The "genetic model" is de- 
scribed in section III. Experimental data are taken from the 
CRC Handbook. 15 £ gap is the HOMO-LUMO gap. 

S ao v\ V2 V3 oh 2 o E coh E sa , p 

model (A' 5 ) force constants (au) (A 3 ) (Ry) (Ry) 

point 0.33 1.029 0.099 1.002 2.3 0.77 0.81 

dipole 0.45 0.25 1.029 0.104 0.840 1.8 0.82 1.03 

genetic 0.47 0.27 1.029 0.065 1.061 1.5 0.76 0.66 

exp. 1.029 0.100 1.062 1.4 0.75 0.82 



FIG. 2. A cartoon of the water dimer. Strictly the two hydro- 
gen atoms at the left are eclipsed and /3 is the angle between 
the O-O bond and the plane of the left hand molecule. In the 
instance that a = we talk of a linear hydrogen bond. 




each model we show the charge transfer, S, larger in 
the dipole model since its derived dipole moment is to 
be partially canceled by the induced moment to achieve 
the target experimental value. The symmetric force con- 
stant associated 16 with the vibrational frequency v\ is 
fitted in the pair potential, but the angular and asym- 
metric frequencies are predictions. First we note that v-i 
is well rendered and we remind the readers that the an- 
gular forces in TB arise from the spa matrix elements, 
transforming like the Slater-Koster table. 18 The predic- 
tions of the asymmetric force constant are revealing. It 
is observed, and also predicted in the local density func- 
tional approximation, that the asymmetric mode is stiffer 
than the symmetric. Intuitively one expects the opposite: 
since the symmetry is broken there should be greater op- 
portunity to relax the electronic structure and hence the 
curvature of the energy-stretch curve ought to be smaller. 
In fact this wrong result is predicted in our models and 
the effect is exacerbated in the dipole model for which the 
direction of anion polarization is free to rotate and hence 
lower the energy. Therefore both models at the intuitive 
level invert the ordering of v\ and v 3 and the dipole model 
is of course the worse culprit. The total polarizability of 
the monomer is well known and can be calculated simply 
by applying an electric field and recording the result- 
ing dipole moment. 2 We find that the polarizability is in 
reasonable agreement with experiment, the dipole model 
being quite accurate in this regard. The cohesive energy 
of the monomer with respect to its atoms is rather well 
reproduced in either model. 



B. Dimer 

The curious structure of the water dimer in vacuo is 
rather well known. We show in Fig. 2 a cartoon of the 
geometry with the principal length and bond angles indi- 
cated. In order to adapt our intuitive models to describe 
the dimer we need to make TB parameters for the 0-0 
interactions. We note that the 0-0 distance, marked 
i?oo m Fig. 2 is not much different from that in a typ- 
ical ceramic oxide, such as zirconia; therefore as a first 
guess, with some modifications we try the hopping pa- 



rameters that we previously used in molecular statics and 
dynamics calculations in that material. 4,19 In that work 
we employed no 0-0 pair potential since these anions 
were second neighbours in the fluorite lattice. For the in- 
tuitive model for water we choose a simple pair potential 
of the form 

tp(r) = Ar~ n e- pr (6) 

having A, n and p all greater than or equal to zero. 
We adjusted these three parameters by eye to obtain 
a reasonable fit of the energy versus 0-0 bond length 
(at fixed angles f3 and a) calculated with the GAMESS 
program 20,21 using coupled clusters at the CCSD(T) 
level. We can then ask, how well do the intutitive mod- 
els predict the angles (3 and a shown in Fig. 2 which 
have been calculated by Klopper et al. 22 and found to 
be a = 5.5° and /3 = 124.4°, while R 00 = 5.50 bohr. 
We use molecular statics to relax the water dimer and 
we find in the point charge model, i?oo = 5-34 bohr, 
a = 13.13° and /3 = 172.7°; while the dipole model 
predicts , R 00 = 5.50 bohr, a = 5.9° and /3 = 102.5° 
This is very encouraging. First the models do not pre- 
dict that the two molecules line up with their monomer 
dipolcs aligned antiparallel — a configuration one might 
well regard as natural from an elementary electrostatic 
viewpoint, and one that DFT in the generalized gradient 
approximation (GGA) predicts to be remarkably close 
(but higher) in energy than the observed orientation of 
Fig. 2. Second, the dipole model at least makes a good 
rendering of the angles /3 and a considering the subtlety 
of the structure. Furthermore the fact that the dipole 
model is evidently superior to the point charge model in 
this regard argues strongly for the correctness in accou- 
ting for the anion polarizability in the construction of a 
TB model for water. 



III. FITTED MODELS 

We hope that the reader is now convinced that the self 
consistent polarizable ion TB theory provides a proper 
framework for the description of the structure, bonding 
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TABLE II. Parameters of our three models — intuitive point 
charge and dipole models and the model fitted by genetic 
evolution. See the text for meanings of the various symbols. 
All quantities are in atomic Rydberg units. 

model: point dipole genetic 



On site parameters 



£s (H) 


-1 


-1 


-1 


£s (O) 


-2 


-1.45 


-1.51 


% (O) 


-1.038 


-1.038 


-1.20 


[7(H) 


1 


1 


1.08 


U(0) 


0.885 


0.77 


1.16 


^^spp 




-0.84 


-0.924 



O-H hopping integrals 





-0.428 


-0.428 


-0.348 


Vspa 


0.550 


0.550 


0.313 




2 


2 


1.48 




2 


2 


1.98 


n c 


4 


4 


4.04 


r c /r 


1.8 


1.8 


1.92 


O-H pair potential 


A 


0.848 


0.937 


0.552 


m 


3.190 


3.045 


3.362 


m c 


6 


6 


6.04 


r c 


2.9 


2.9 


3.04 



O-O hopping integrals 



Vssa 


-0.072 


-0.072 


-0.080 


Vspa 


0.084 


0.084 


0.050 


Vppa 


0.06 


0.06 


0.00012 




-0.01 


-0.01 


-0.004 


n S ser 


2 


2 


2 


l^spcr 


2 


2 


2 




3 


3 


3 




3 


3 


3 


n c 


6 


6 


4 


d Q 


5.6 


5.6 


5 


r c 


9 


9 


6.8 



O-O pair potential 



A 


10 5 


1.5 x 10 s 




n 


9.7 


6 




P 





1.2 










0.010 


u 2 






0.647 


ro 






5.992 


r\ 






5.494 


r c 






6.110 



TABLE III. Properties predicted, or fitted, of the water dimer. 
Target values are taken from Klopper et al. 22 

model: point dipole genetic target 

Roo (bohr) 5.3423 5.5011 5.5091 5.4991 

a (deg.) 13.1 5.9 3.0 5.5 

/3 (deg.) 172.7 102.5 113.7 124.4 

E coh (mRy) -18.2 -16.8 -15.1 -15.9 



FIG. 3. Cohesive energy verus O-O distance in the dimer at 
fixed a = 5.5° and /3 = 124.4°, using the two intuitive and the 
fitted models, compared to quantum chemical calculations at 
the CCSD(T) level using the GAMESS suite of programs. 20 ' 21 




4 5 6 7 

d 00 (bohr) 



and energetics of the water monomer and dimer. In that 
case there should be no objection to the next develop- 
ment which is to seek refined hamiltonian matrix ele- 
ments, pair potentials and their bond length dependences 
via a thorough search in the parameter space. To this 
end, we have adopted Schwefel's multimembered evolu- 
tion strategy. 23,24 As we have seen in section II a set of 
TB parameters may be used to predict certain proper- 
ties. Indeed any set of parameters will give rise to values 
of a number of chosen properties which can be compared 
to target values taken from experiment or coupled cluster 
calculations. A suitably weighted sum of squares of dif- 
ferences becomes an objective function and the aim is to 
find that set of parameters which gives the smallest value 
of the objective function, subject to certain contraints 
on the parameters so that they continue to be physically 
motivated. This last point is very important — we cannot 
accept a model whose parameters violate certain basic 
truths; for example e p — e s > 0, V ssa < and so on. 
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We also have an instinct for the approximate sizes and 
appropriate scaling laws for the hopping integrals which 
are informed by canonical band theory. 9 The procedure 
is this. A computer script is written that, for a given set 
of TB parameters, calculates the required set of proper- 
ties and evaluates the objective function. The evolution 
algorithm repeatedly calls this script with TB parame- 
ters that belong to successive parents and offspring until 
a global minimum is hopefully found. In the process, a 
huge set of parameters is tested and we are free to browse 
this set to select attractive looking models and to inves- 
tigate how certain properties are governed by particular 
parameters. 



A. Monomer 

In the case of the monomer, the bond length, dipolc 
moment, frequencies and polarizability are used in the 
fitting and the resulting model is added here to table I. 
We note that the correct ordering of v\ and v% is acheived, 
at the expense of an angular force constant that is sig- 
nificantly too small. We believe that for the simulations 
of liquid water, to be described in section VI below, the 
polarizability of the monomer plays a large role and it is 
therefore important at this stage that this is now prop- 
erly rendered in the "genetic model" , so called. 



B. Dimer 

It is a feature of our "ground up" philosophy towards 
the construction of TB models that, having established 
a model for the monomer, we do not further adjust its 
parameters; and continue to the fitting of the oxygen- 
oxygen interactions only to the properties of the dimer. 
Our objective function is constructed using the relaxed 
geometry of the dimer, in particular, RoOi & an d ft ob- 
tained from a full geometry optimization. In addition we 
do calculations of the cohesive energy at a set of 0-0 
distances at fixed values of a and /?, namely their target 
values as found in the literature. 22 Five such datapoints 
are included in the objective function. Fig. 3 shows the 
energy that we have calculated 20 ' 21 in coupled clusters at 
the CCSD(T) level compared to data points from the two 
intuitive models and our final "genetic" model. We note 
that the exponential pair potential, (6), which is typical 
for standard solid state TB models is too steeply repul- 
sive compared to the CCSD(T) curve. This is clearly a 
peculiar signature of the hydrogen bond. We have chosen 
a weaker pair potential based on the very first, quadratic, 
pair potential proposed by Jim Chadi to describe the har- 
monic properties of semiconductor s — p bonding. 25 

ip(r) = C/ ie + [/ 2 e 2 

where e = (r — r a ) /r is the fractional change in bond 
length relative to some reference distance, r$. Of course a 



parabola does not describe the CCSD(T) curve in Fig. 3 
at large 0-0 distances; we replace the potential <p(r) 
with a fifth degree polynomial beyond some distance, n , 
which is constructed to be continuous and twice diffcr- 
cntiable at n and to vanish smoothly at a cut-off dis- 
tance, r c . This is consistent with modern practice in the 
bond order potentials. 26 The behavior of this curve at 
long distances is dictated by electrostatics and van der 
Waals interactions. The latter are explicitly absent in 
our model while the electrostatics are included through 
the self consistent hamiltonian (section II A 4). 

All of the parameters of our three models are gathered 
into table II. The properties of the dimer are displayed 
in table III. 



C. The hydrogen bond 

It is worthwhile to make a few preliminary remarks 
concerning the hydrogen bond at this stage of the devel- 
opment. In a sense the water dimer in its geometry of 
Fig. 2 is an archetypal example; and one may ask what is 
the nature of such a bond and what determines the values 
of i?oo, oi and /?? The repulsive force which we model 
with a pair potential must arise in real life from the closed 
shell repulsion and valence-core overlap. Our view is that 
the principal attractive force comes about from the mu- 
tual static polarizability of the two monomers and that 
in addition the point charge electrostatics plays a role in 
the angular disposition of the molecules. On the other 
hand we recall again that the 0-0 distance is only mod- 
erately smaller than that in a typical metal oxide which 
can be modeled as mixed ionic-covalent bonding. 1 For 
this reason we include the sp and pp bonding integrals in 
our model. If we leave these out, we can make a plausible 
model, but we find we cannot reproduce the angles well. 
One can imagine that negotiation between the a and 7r 
bonds is what results in the final choice of the angle (3. 



IV. THE HEXAMER 

The nature of the hydrogen bond as manifested in the 
dimer can be further studied in the structure and prop- 
erties of the hexamer in vacuo. Four isomeric structures 
have been identified, 27 and remarkably these are practi- 
cally degenerate in energy even though their shapes and 
the number of hydrogen bonds are very different. 27 The 
four isomers are displayed in Fig. 4 which shows the struc- 
tures after geometry optimzation using the genetic fitted 
TB model of section III. As in the dimer, the hydrogen 
bonds are not linear: again one sees the angle equiv- 
alent to a in the dimer is some 5°. The exception to 
this is the cyclic hexamer which shows almost exactly 
linear hydrogen bonds. Significantly this is the least sta- 
ble of the four. Although the TB model overestimates 
the energy differences with respect to the most stable, 
the prism, by as much as an order of magnitude it is 



7 



FIG. 4. Structures of the six hexamer isomers considered in 
the text, relaxed using the TB model. The labeling of the 
oxygen atoms is such that the O-O bond lengths can be read 
off in table V 




very significant that the ordering predicted using quan- 
tum chemistry at the MP2 level 27 is reproduced by our 
model. This is especially noteworthy since the GGA pre- 
dicts the ordering quite reversed, namely in that approxi- 
mation, the book is most stable followed by the cyclic, the 
cage and prism in that order. 27 Our TB model predicts 
that the binding energy of the prism hexamer, relative 
to the isolated monomers is —0.132 Ry. On the other 
hand the binding energy relative to three isolated dimers 
is —0.087 Ry. This reflects the additional binding due 
first to the establishment of three more hydrogen bonds; 
and second, in our view having greater significance, the 
attraction due to the mutual polarization of the water 
monomers in forming the hexamer from dimers. In fact 
we have found, in common with previous first principles 
calculations, 28 ' 29 that the dipole moments of component 
monomers increase as they are assembled into succes- 
sively larger clusters. This is shown in table IV. We will 
return to this point when we discuss the structure of liq- 
uid water. We should note that unlike density functional 
theory, the tight binding model provides an unambigu- 
ous measure of individual molecular multipole moments; 
this does not however detract from the fact that these 
are, in principle, unmeasurable and therefore ill-defined 
quantities. 

The relaxed structures are compared with the MP2 cal- 
culations by displaying the 0-0 hydrogen bond lengths 
in table V. Our relaxed structures compare extremely 
well with the quantum chemistry results; and we remark 
that the good agreement and the improvement in our 
model over the GGA is because we have rather care- 



TABLE IV. Dipole moments on individual water molecules 
found in the monomer, dimer, trimer . . . hexamer using our 
fitted genetic, and intuitive point charge TB models and com- 
pared to available similar results using the GGA to density 
functional theory. Note in connection with remarks in sec- 
tion VI that the point charge model does display the expected 
increase in dipole moment with size of cluster, but this is non 
monotonic, the trimer being an exception to the trend, and 
the magnitude of the effect is probably exaggerated by the 
anomalously large polarizability of the monomer (Table I). 
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1.86 


2.07 


2.14 


2.21 


2.31 


2.44 • • 


• 3.28 


GGA 


1.87 


2.1 


2.4 








• 2.95 



fully obtained agreement over a wide range of 0-0 bond 
lengths with the CCSD(T) curve in Fig. 3. It is well 
known that the local density functional theory, even in 
gradient corrected form, has difficulty with the energet- 
ics of the hydrogen bond and so we have been able to 
leap-frog these problems in our TB model by comparing 
to the CCSD(T) not the GGA. 



V. THE STRUCTURE OF ICE 

We now turn to the solid state to examine how our 
model reproduces the structure of ice. This is not the 
place to make an exhaustive study of the many solid state 
phases, nor to attempt a prediction of the phase diagram, 
we will leave that to future work. Most solid phases have 
disordered hydrogen atom sites, but the phase Ice- XI is 
thought to be "proton-ordered" so we test our model 
against this structure. Hirsch and Ojamae 30 have iden- 
tified a number of putative structures of Ice Ih of which 
Ice- XI is found to have lowest energy. In Fig. 5 we show 
the structure of Ice- XI after geometry optimization us- 
ing our genetic TB model. We find its cohesive energy 
with respect to the isolated water monomers is 39 mRy 
per monomer, or 12.2 kcal/mol. This result is consistent 
with the experimental lattice energy of 13.3 kcal/mol. 30 
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TABLE V. Properties of the six hexamer isomers considered 
in the text after geometry optimization using the TB model. 
Energy differences, AE, are with respect to the Prism isomer. 
The results labeled "MP2" are taken from Santra et al. 27 



Prism Cage Book Cyclic 

TB MP 2 TB MP 2 TB MP2 TB MP2 

AE (mRy) 2.8 0.2 6.9 1.5 11.8 5.0 

0- O bond lengths (bohr): 

1- 2 4.975 5.280 5.055 5.205 4.970 5.094 5.009 5.134 

2- 3 5.639 5.522 4.914 5.047 4.955 5.093 5.009 5.134 

3- 1 5.597 5.545 

3- 4 5.208 5.265 5.045 5.298 5.009 5.134 

4- 5 5.382 5.240 5.003 5.225 5.061 5.249 5.008 5.134 

5- 6 5.564 5.424 5.542 5.520 5.299 5.261 5.009 5.134 

6- 4 5.627 5.556 

1- 4 5.533 5.512 

2- 5 4.926 5.028 5.570 5.573 

3- 6 5.013 5.185 5.566 5.494 5.520 5.544 

6-1 4.998 5.149 4.972 5.081 5.009 5.134 



TABLE VI. Lattice constants and density of ice-XI at 0°K 
calculated in our genetic and intuitive point charge TB models 
compared to GGA predictions and experiment. Densities of 
liquid water are also shown for comparison (in g/cm 3 ), and 
discussed below in section VI. BLYP 32 ' 33 and PBE 34 denote 
specific gradient corrected density functionals. 







Ice 


-XI 




Liquid 




a 


b 


c 


P 


P 


TB (genetic) 


8.422 


14.506 


13.668 


0.967 


0.926 


TB (point) 


8.488 


14.754 


13.900 


0.928 


1.030 


BLYP 


8.227 a 


14.394 a 


13.619 a 


0.995 a 


0.80 b 


PBE 


8.282 c 


14.405 c 


13.536 c 


1.000 c 


0.88 b 


experiment 


8.438 d 


14.850 d 


13.780 d 


0.935 d 


1.00 



a Reference [30 . 
b Reference [35 . 
c Reference [36 . 
d Measurements at 5°K. 



We also find that the optimized structure of the proton- 
ordered Ice-Ih "phase number 6" of Hirsch and Ojamae 30 
has an energy higher than Ice-XI by 0.16 kcal/mol which 
may be compared with the DMol3/BLYP calculation re- 
sult of 0.05 kcal/mol. 30 It is significant that classical force 
field models find, conversely, that "structure number 6" 
is more stable than Ice-XI. Our optimized Ice-XI has a 
density at 0°K of 0.97 g/cm 3 compared to the experimen- 
tal value of 0.93 g/cm 3 at 5°K. 31 On the other hand the 
density of Ice-XI as calculated in the GGA is as large as 
1.00 g/cm 3 as we report in table VI which also shows our 
calculated lattice constants of ice-XI compared to GGA 
calculations and experiment. 



VI. THE STRUCTURE AND PROPERTIES OF LIQUID 
WATER 

Finally we apply our TB models to predict properties 
of liquid water. For this we take a cubic box of side 
29.75 bohr containing 128 water molecules, having pe- 
riodic boundary conditions and previously equilibrated 
in molecular dynamics (MD) using the SIESTA program 
and the BLYP GGA. 44 We further equilibrated for 100 ps 
in an NVT ensemble at 27 C and 1 g/cm 3 before an ad- 
ditional 100 ps to gather statistics in an NVE ensemble. 
We have also made a simulation in an NPT ensemble at 
27 C and zero external pressure for a total of 100 ps, tak- 
ing statistics from the last 50 ps of simulation time. We 
employ the reversible integrators with Liouville operators 
developed by Martyna et al. 45 using just a single Nose- 



TABLE VII. Some properties of liquid water calculated by our 
intuitive point charge and genetic dipole models, compared 
to experiment and published GGA calculations. p to t and p ln d 
are total and induced average dipole moments, taken from 
Fig. 7, e(0) is the static dielectric constant and D ae u is the 
self diffusion coefficient. PBE and BLYP have the meanings 
of table VI. 





Ptot Pind 


e(0) 


Aself 


model 


(D) (D) 




10~ 5 cm 2 /s 


experiment 


2.95±0.2 a 


78, b 68 c 


2.2 d 


genetic (27° C) 


3.28 0.74 


86.7 


3.0 


point (27° C) 


2.51 


58.7 


3.5 


PBE 


2.95, c 3.09 f 


67±6, f 75 s 


1.6 f 


BLYP 






0.25, h 0.55 1 



a 27° C, Reference [38]. 

b 27° C, Reference [39 and 40]. 

c 57° C, Reference [40]. 

d 27° C, Reference [41]. 

° 45° C, Reference [29]. 

f 57° C, Reference [42]. 

g Extraplolatcd to 27 C from data in reference [42]. 
h 35° C, Reference [43]. 
1 32° C, Reference [44]. 
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FIG. 5. The structure of ice-XI after geometry optimization 
within our fitted genetic TB model. A unit cell is marked by 
straight lines. 




Hoover chain and particle and barostat relaxation times 
of 100 fs and 1 ps respectively. We use atomic masses of 
oxygen and hydrogen (not deuterium) and time steps of 
1 fs in NVT and NVE and 0.25 fs in NPT ensembles. 

First we show radial distribution functions (RDF) in 
Fig. 6. These show the same sorts of features in compari- 
son to experiment as found in previous calculations using 
GGA. 35,43 In particular the first neighbor O-H peak is 
narrower than experiment, due probably to the quantum 
nature of protons that is usually neglected in MD simula- 
tions. The most critical RDF is that of 0-0 bonds since 
this describes the complex of hydrogen bonds responsi- 
ble for many of the subtle properties of water. The RDF 
from our intuitive point charge model is very similar to a 
recent calculation 46 using the the SCC-DFTB method, 7 
which as discussed in section II above is equivalent in con- 
struction to our point charge model. The point charge 
model suffers from a nearest neighbor 0-0 peak which 
is too sharp and too high and a lack of structure in the 
RDF beyond first neighbors; that is to say, overstruc- 
tured in the first solvation shell and understructured in 
more distance shells. Both these shortcomings are reme- 
died by our genetic dipole model at the expense of the 
first neighbour 0-0 peak falling beyond the experiment 
indicating an average hydrogen bond length that is too 
long. It is notable that this average length is quite in- 
sensitive to the simulation box volume, since the NVE is 
done at a density of 1 g/cm 3 , while the average density 
in the NPT simulation was 0.926 g/cm 3 (see table VI). 

We use Figs. 7 and 8 to illustrate the point that we 
have made earlier (see Tables IV and VII) that the molec- 
ular dipole moment increases significantly as the water is 
built up in clusters from individual molecules having a 
dipole moment of 1.86 D to a total dipole moment in liq- 
uid water of about 3 D. It is a particular feature of our 
dipole model that this effect is properly captured, in con- 
strast to our point charge model in which the molecular 
dipole moments do not grow monotonically with clus- 
ter size (see Table IV) and do not reach the expected 



FIG. 6. Calculated radial distribution functions in liquid 
water, compared to experimental neutron diffraction data. 37 
NVE and NPT refer to our genetic dipole model simulations 
in microcanonical and canonical ensembles. NVE is done at 
a density of 1 g/cm 3 . The average density in NPT simulation 
is displayed in table VI. We also show NVE simulations in 
our intuitive point charge model. The lowest panel shows the 
integrated quantity m(r) = J Q r goo(r')r' 2 dr' and we mark the 
position of the first minimum to show the average number of 
0-0 bonds. As is generally known this is larger than the 
number, 4, associated with the hexagonal ices, such as ice-XI 
shown in Fig. 5. 
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FIG. 7. Distribution of the magnitude of the dipole monent 
in an NVE molecular dynamics. In the case of the intuitive 
point charge model there is just the dipole moment arising 
from point charge transfer, indicated by a dotted line. In the 
genetic dipole model the point charge contribution is larger, 
consistent with the picture for the monomer, shown in Fig. 1; 
however, unlike in the monomer, the partial cancelation from 
the oppositely oriented induced dipole moment is not com- 
plete (see Fig. 8) so that the total dipole moment is on av- 
erage larger than in the point charge model and is consistent 
with measurement and GGA calculations. 
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~3 D in the bulk limit. Moreover most of the increase in 
dipole moment is probably a result of the polarizability 
of the monomer being nearly a factor of two too large 
in the intuitive point charge model (see Table I). The 
same failing, and to a much greater degree, is apparent 
in semiempirical schemes such as AMI and PM3 (see ta- 
ble 4 of ref [47] ) . The quantitative rendering of the dipole 
moments as a function of environment as evidenced in 
Table IV is probably the one principal advantage con- 
veyed by using a model with polarizable ions. Classical 
models usually establish a molecular dipole moment of 
the monomer to about 2-3 D at the outset of the con- 
struction of the model. 

Table VII lists some properties of liquid water arising 
from our simulations, compared to previous GGA cal- 
culations and experimental results. The static dielectric 
constant is calculated using the method of Sharma et 
al. 42 and the self diffusion coefficient is extracted as the 
linear slope of the mean square displacement. We should 
note that both models give a reasonable prediction of 
the dielectric constant, the genetic being rather better, 
which we attribute to the anion polarizability. The intu- 
itive point charge model overestimates D se if which is con- 
sistent with recent SCC-DFTB simulations, 46 although 
these latter overestimate by a great deal more, finding 
D se if = 11.1 ± 0.4 and 6.5 ± 0.2 x 10~ 5 cm 2 /s in two 
different modifications of the SCC-DFTB model, which 
are very much greater than our 3.5 x 10~ 5 cm 2 /s. This 



FIG. 8. A cartoon to illustrate the incomplete cancelation 
of the point charge transfer dipole moment in liquid water. 
The left hand sketch shows molecules in the point charge 
model whose dipoles are necessarily aligned along the cen- 
terline of the molecule and whose magnitude is close to that 
in the monomer. In the dipole model with induced dipoles, 
the right hand cartoon explains how the molecular dipole mo- 
ment is increased: in the monomer (see Fig. 1) the induced 
dipole moment is constrained to point opposite to the point 
charge transfer moment, leading to a reduction in the total 
moment. Conversely in the liquid (and also in the solid and 
clusters discussed above) the induced moment "measures" the 
total electric field at the anion, not just that due to the two 
closest hydrogen ions; in this way the partial cancelation of 
the monomer is supressed and the resulting total moment is 
closer to that of the point charge transfer which is larger than 
the corresponding quantity in the point charge model. 
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is very encouraging, firstly indicating that our intuitive 
point charge model is better in this regard than the SCC- 
DFTB model; more importantly by extending the theory 
to allow anion polarizability results in a significant im- 
provement in the rendering of the self diffusion in the 
liquid. In contrast to the GGA method, 42 we reproduce 
D sc \[ quite accurately without recourse to artificial eleva- 
tion of the temperature. 



VII. DISCUSSION 

A. General remarks on the new model 

The central result of this work is that it is possible 
to construct TB models of water in a "ground up" ap- 
proach; that is, by fitting parameters to observed and 
calculated properties of the monomer and dimer alone 
the model is found to be sufficiently robust and trans- 
ferable to be able to predict many of the known prop- 
erties of water. We believe that the rather transparent 
construction that we have outlined in section II allows 
the models to provide some insight into the very subtle 
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properties of water, many of these being manifestations 
of the hydrogen bond. In section III C we made some 
observations on the nature of the hydrogen bond, assert- 
ing that its principal two components are an attractive 
force arising from the mutual static polarizability of the 
monomers countered by a short ranged repulsion which 
we model with the usual pair potential of TB theory. Of 
secondary importance is the additional attraction due to 
orbital overlap which is modeled by hamiltonian hopping 
matrix elements as in all electronic structure localized 
orbitals methods, including the density function descrip- 
tion. We have not addressed the question of van der 
Waals forces. While we expect these to be weaker and 
to decay more rapidly with 0-0 distance than the static 
polarizability, there is recent evidence from DFT stud- 
ies that their inclusion significantly improves upon the 
local density approximation. 27,35 One could argue that 
by careful fitting to the CCSD(T) curve in Fig. 3 we 
have "renormalizcd" our pair potential to include some 
"l/r 6 -like" attraction as would be done in a classical po- 
tential. Conversely, as in the DFT, we could extend our 
TB theory to include dynamic polarizability explicitly. 14 
Whether this can lead to an improved description of the 
hydrogen bond must be left for future work. 



B. The density of water and ice 

The most striking and anomalous feature of water is 
the fact that ice is less dense than its liquid, leading to the 
possibility of life on earth. Therefore we would very much 
like our models to display this. Unfortunately as we have 
seen in table VI our genetic dipole model gives an ice den- 
sity of 0.97 g/cm 3 at 0°K. Wc find in an NPT simulation 
at -3 C that this decreases somewhat to 0.96 g/cm 3 but 
this is still greater than the average density of the liquid 
in the NPT simulations at 27 C, namely 0.926 g/cm 3 . 
On the other hand this situation is nothing like as bad 
as the GGA densities reported in the literature. 30,35 ' 36 
These find the density of ice-XI to be as much as 1 g/cm 3 
while the liquid density is as low as 0.8 g/cm 3 in the 
BLYP functional and 0.88 g/cm 3 in the PBE functional. 
This situation is rescued using a functional including dy- 
namic van der Waals forces but at the expense of a poor 
rendering of the 0-0 RDF. 35 From the DFT point of 
view the question of the densities of ice and liquid water 
is clearly not yet resolved. In the TB framework it is 
most surprising and at the same time very encouraging 
that the intuitive point charge model correctly predicts 
that ice will float on water. This raises two questions. 
(i) Is the inclusion of anion polarization an unnecessary 
and uncontrollable complexity? In answer, we can cer- 
tainly see from the above that the genetic dipole model 
performs better than the point charge model in respect 
of the RDF and properties of the monomer and dimcr. 
On the other hand we have not refined the point charge 
model but have kept it at the "intuitive" level, (ii) Could 
a point charge model be constructed that will adopt the 



TABLE VIII. O-O bond lengths in the structures studied here 
using our two models (in atomic units of length). We have 
used the cyclic hexamer as an example. In the case of ice-XI 
and liquid water these are averages. 

model dimer hexamer ice-XI liquid 
genetic 5.51 5.00 5.13 5.52 
point 5.34 5.01 5.22 5.26 



benefits of the dipole model and retain the correct differ- 
ence between the liquid and ice densities? This must be 
a subject of further research but we suspect not. First we 
know from the failings of the SCC-DFTB simulations, 46 
in respect of both 0-0 RDF and diffusivity that the TB 
description can fail at the point charge level, although 
our intuitive point charge model corrects the worse short- 
comings of the SCC-DFTB. Second, and this is a point 
we wish to reemphasize, it is only by including anion po- 
larizability that the trend of increasing molecular dipole 
moment with cluster size into the limit of the solid and 
liquid can be reproduced (table IV). We believe that it is 
this one element that elevates our TB model into direct 
competition with the very much more (computer) time 
consuming DFT models. 

We can use some of the results above to speculate upon 
the question, why is ice less dense than water? The con- 
ventional wisdom has it that hexagonal ice having four 
hydrogen bonds to each molecule is less densely packed 
than the liquid in which, as we see from the lowest panel 
in Fig. 6, each molecule has on average between five and 
six neighbouring molecules. However this cannot be the 
whole story because this is true also in our genetic model 
and indeed published DFT simulations in which ice is 
(wrongly) more dense than water. It is important to 
note that, in contrast to the solid state, the average 0-0 
bond length in liquid water in our models depends on 
the model parameters, but not strongly on the size of 
the confining box, that is, on the volume. This is evident 
from the first peaks of goo{ r ) from the NVE and NPT 
simulations in Fig. 6. Some light is thrown on the mat- 
ter by table VIII which shows (average) hydrogen bond 
lengths from our calculations using the genetic dipole and 
intuitive point charge models. The point charge model 
gives the bond length and hence density of ice-XI much 
more accurately compared to experiment than the ge- 
netic dipole model. Furthermore since the first peak of 
9oo{ r ) hi the intuitive point charge model is close to the 
experiment (Fig. 6) the consequence is that this model, 
perhaps fortuitously, predicts the right ordering of the 
densities of ice and water. This may well indicate that 
this will be correctly rendered by any model that simul- 
taneously predicts both the density of ice and the position 
of the first peak of goo(r)- However it is not obvious that 
such a model can be easily found without specifically fit- 
ting these two properties — this is of course possible for 
classical models. 
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VIII. CONCLUSIONS 



are the Gaunt integrals for real spherical harmonics. 14 



We have demonstrated models for water at "intuiti- 
tive" and fitted levels in the point charge transfer and 
dipole charge tight binding approximations. They have 
the benefit of quite clear physical and chemical insight by 
virtue of the way they are constructed, and at the same 
time possess quite remarkable predictive power. This is 
evident from the "ground up" approach that we have 
adopted in which only properties of the monomer (in re- 
spect of O-H interactions) and the dimcr (in the model- 
ing of the hydrogen bond) are fitted. Properties of clus- 
ters, ice and liquid water are then predicted with good 
quantitative agreement with experiment and published 
DFT calculations. Indeed in some notable respects our 
models are superior to the DFT, in particular in the cohe- 
sive energy curve of the dimer (Fig. 3) and in the density 
of ice-XI and the liquid (table VI) which are particularly 
poorly described by the GGA. 

We expect our model to find a niche where nanosecond 
simulations of some thousands of particles are called for 
while retaining a proper self consistent quantum mechan- 
ical description of the chemical bond. In particular, we 
believe that the model will provide insight into the com- 
plexities of the hydrogen bond and in the role of polar 
protic solvents in chemical reactions. 
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Appendix A: Generalized Madelung matrix 

The generalized Madelung matrix B is proportional to 
the KKR structure constant matrix, B 48 > 49 



Br'l' rl — 



47T 



{21 + l)!!(2f + 1)!! 



Bt)i 



R'L' RL 



where 



Br<l<rl = 4tt^(-1)^" - 1)!! Cvlv> K L »(R, — R') 



and 

K L (r)=r- e - 1 Y L (r) 
is the solid Hankcl function. 

Cwtjl = J J df2 Yjjt Ytj Yi 
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